Necesitamos un editor

Más práctico: ejecutar múltiples operaciones simultáneamente.
Para ello: se escriben todas en un fichero de texto simple, pero con extensión .r en lugar de .txt.
Es suficiente un editor de texto simple. En Windows: Bloc de notas (en OS X: TextEdit). Los editores con formatos complejos como Microsoft Word no valen.
Otros editores facilitan la introducción y ejecución de las instrucciones.
Ej.: Notepad++, Sublime Text, Visual Studio Code, Geany, Nano, Emacs, VIM, …
Enfocado al uso de R y el más usado para ello: RStudio (nombre del editor y de la empresa que lo produce y comercializa).
- RStudio es una interfaz que tiene el objetivo de facilitar y agilizar el uso de R y otros formatos que interactuan con R (documentos pdf, archivos word, páginas y aplicaciones web, etc., todos ellos con código R ejecutado).
- En particular ofrece menús para especificar el directorio de trabajo, instalar, actualizar y/o cargar librerías, y abrir y visualizar bases de datos, entre otros. Asímismo dispone de atajos de teclado para introducir los comandos más habituales de manera más rápida.
- RStudio no es software libre, aunque su edición básica es de código abierto y gratuita.
- https://www.rstudio.com/

Permite tener un fichero con extensión .r en el que se puede escribir el código (con autocompletado) y ejecutarlo
Permite añadir instrucciones mediante atajos de teclado

Algunas nociones sobre RStudio

Las instrucciones en el R-script en RStudio se ejecutan con Control+Enter

Qué es lo primero que nos encontramos en R? Objetos, funciones, comentarios, (directorios) y librerías

Objetos elementales:
- Numéricos, como en una calculadora (vimos algunos hace unas diapositivas)
- Cadenas de texto, siempre limitadas por comillas simples o dobles, desde un carácter, por ejemplo "a", o 'z', hasta un texto largo como "Lorem ipsum dolor sit amet..." o 'Lorem ipsum dolor sit amet...'
Funciones: ejecutan las instrucciones que las definen. Genéricamente se ejecutan escribiendo el nombre de una función f seguido de los objetos o valores que necesita entre paréntesis.
- Cuando no es necesario ningún valor, estos paréntesis se escriben sin incluír nada, así: f()
- En la mayoría de casos, necesitan evaluar objetos. La evaluación de un objeto x con una función f se codifica f(x).
Un comentario es un fragmento de código, texto… no ejecutable. Se señalan situando a la izquierda un símbolo #
Directorio de trabajo: el entorno de R tiene como referencia un directorio concreto. Esto nos permite no tener que introducir la dirección de la carpeta cuando queremos acceder a algún elemento de esta (por ejemplo, ficheros con bases de datos).
- Para saber cuál es el directorio de trabajo actual utilizamos la función getwd;
- Para cambiar el directorio se utiliza la función setwd evaluando la dirección del directorio de trabajo en que queremos situarnos ahora;
- Se puede saber qué archivos hay en el directorio donde R está trabajando con la función list.files;

getwd() # el directorio donde R está trabajando
setwd("P:/projects/Curso-R") # (en Windows) se tienen que sustituir los \ por /
list.files()

O bien, mediante RStudio:

`Session > Set Working Directory > Choose Directory`
`Session > Set Working Directory > To Source File Location`

Librerías:
- Para instalar una librería llamada tidyverse evaluamos con la función install.packages la cadena consistente en el nombre de la librería, "tidyverse";
- Para cargar una librería llamada readxl, evaluamos con la función library (o bien require, una función distinta), bien la cadena con el nombre de la librería "readxl" o bien (el nombre de) la librería misma, readxl.

install.packages("tidyverse") # instala la librería tidyverse, que a su vez instala las librerías readxl, dplyr, tidyr, haven, y otras (https://www.tidyverse.org, https://github.com/tidyverse/tidyverse)
library(readxl)
# library("readxl")
# require(readxl)
# require("readxl")
# (.packages()) # vemos qué librerías tenemos cargadas en este momento

O bien, mediante RStudio:

Panel `Packages > Install`

Siguientes pasos: formas de obtener ayuda, argumentos de una función y abrir un fichero con formato de Excel

R dispone de una ayuda muy completa sobre todas las funciones, procedimientos y elementos que configuran el lenguaje
Se puede acceder a ella con la función help() o mediante ?. Los siguientes ejemplos nos muestran los documentos de ayuda de las respectivas funciones.

help(install.packages)
?library

Para evaluar un objeto mediante una función, implícitamente estamos asignando el objeto a un argumento de la función. Por ejemplo, al ejecutar library(readxl) asignamos (el nombre de) la librería readxl al primer argumento de la función library, que es package. Entonces, otra forma de ejecutar library(readxl) es ejecutar library(package = readxl). En este caso, es exactamente lo mismo, pero cuando queremos usar otro argumento de una función distinto del primero en general tenemos que nombrarlo.
Además de las funciones, entre otros objetos de los que podemos ver documentos de ayuda están las librerías.

?help
help(package = "readxl") # Informémonos acerca de la librería readxl
?read_excel # ayuda de la función read_excel

Para abrir ficheros de excel disponemos de diversas funciones en las librerías readxl y openxlsx entre muchas otras. Aunque la segunda librería tiene más funciones y opciones para manipular y escribir ficheros .xlsx, la primera dispone de funciones que nos permiten leer también ficheros .xls. En cualquier caso, estas funciones evalúan, bien la dirección absoluta del fichero que hay que abrir, bien la dirección relativa con respecto al directorio de trabajo de R como en el siguiente ejemplo.

read_excel("data_curs_stat/EP1.xls")

## # A tibble: 4,753 × 12
##    q0002_hhid sex   age   mar_1 edu1      phys_hea1                    hea1 dep1  score_lon1 score_sup1 income1 income_inf1
##         <dbl> <chr> <chr> <chr> <chr>     <chr>                       <dbl> <chr>      <dbl>      <dbl>   <dbl> <chr>      
##  1          1 fem   35-49 Yes   Tertiary  None                         70.5 No             3         12      15 Normal     
##  2          2 masc  35-49 Yes   Tertiary  None                         78.9 No             3         12      17 Normal     
##  3          3 masc  65-79 Yes   Primary   None                         66.6 No             3         12      15 Normal     
##  4          4 masc  50-64 Yes   Primary   1 physical health problem    77.4 No             3         11      18 Normal     
##  5          5 fem   50-64 No    Tertiary  None                         52.5 No             3         12      29 Normal     
##  6          6 fem   50-64 No    Tertiary  None                         53.7 No             3         13      25 Normal     
##  7          7 masc  50-64 Yes   Secondary None                         59.1 No             6         10      12 Good       
##  8          8 masc  50-64 Yes   Primary   None                         76.4 No             3         12      14 Normal     
##  9          9 fem   50-64 No    Secondary None                         73.7 No             3         11      26 Good       
## 10         10 fem   80+   No    Primary   2+ physical health problems  28.0 No             6         13      10 Normal     
## # … with 4,743 more rows

Nota: Para abrir ficheros *.csv disponemos de la función read.csv en R, de la función read_csv en la librería readr y de la función fread en la librería data.table entre otras. La librería haven dispone de diversas funciones que permiten abrir y escribir ficheros en formatos de Stata, SPSS y SAS (pueden consultarse en su ayuda, package(help = "haven")). Por otra parte la librería readspss (https://github.com/JanMarvin/readspss) tiene funciones que permiten abrir bases de datos en formato de SPSS encriptadas con contraseña. Podemos encontrar más información en esta breve introducción a R

Abrir bases de datos y ejecutar funciones: R no es como Stata.

La función read_excel abre la base de datos, pero a diferencia de instrucciones como use o import excel en Stata, no la copia a memoria. Así, si en Stata podríamos ejecutar describe o summarize, si en R ejecutamos una función como summary el resultado es la propia función.

summary

## function (object, ...) 
## UseMethod("summary")
## <bytecode: 0x55759c4cd800>
## <environment: namespace:base>

Tenemos que evaluar la base de datos que abrimos con read_excel con la función summary:

summary(read_excel("data_curs_stat/EP1.xls"))

##    q0002_hhid        sex                age               mar_1               edu1            phys_hea1              hea1           dep1             score_lon1     score_sup1       income1       income_inf1       
##  Min.   :    1   Length:4753        Length:4753        Length:4753        Length:4753        Length:4753        Min.   : 0.00   Length:4753        Min.   :3.00   Min.   : 3.00   Min.   : 1.000   Length:4753       
##  1st Qu.: 1327   Class :character   Class :character   Class :character   Class :character   Class :character   1st Qu.:42.83   Class :character   1st Qu.:3.00   1st Qu.:10.00   1st Qu.: 2.000   Class :character  
##  Median : 2615   Mode  :character   Mode  :character   Mode  :character   Mode  :character   Mode  :character   Median :55.64   Mode  :character   Median :3.00   Median :12.00   Median : 5.000   Mode  :character  
##  Mean   : 3145                                                                                                  Mean   :53.97                      Mean   :3.74   Mean   :11.53   Mean   : 8.989                     
##  3rd Qu.: 3807                                                                                                  3rd Qu.:66.06                      3rd Qu.:4.00   3rd Qu.:13.00   3rd Qu.:15.000                     
##  Max.   :92755                                                                                                  Max.   :92.82                      Max.   :9.00   Max.   :14.00   Max.   :35.000                     
##                                                                                                                 NA's   :187                        NA's   :285    NA's   :414     NA's   :1067

Abrir la base de datos mediante la función read_excel cada vez que queremos operar con ella no es práctico. En lugar de esto, guardamos la base de datos en una variable.

Objetos complejos. Clases. Bases de datos

R es un lenguaje de programación que nos permite acceder a objetos con distintas estructuras y manipularlos. Los tipos de objetos más importantes son:

El vector, una colección de datos elementales del mismo tipo básico ordenados.
- Los tipos básicos más comunes son:
  - Numérico
  - Cadena de texto o caracteres
  - Lógico: TRUE, FALSE y NA (missing)
  - Entero, para números exclusivamente enteros
- Habitualmente, los vectores se crean y se unen con la función c
La lista, una colección de objetos de cualquier tipo. Habitualmente se crean mediante la función list
La función
El símbolo, el nombre de un objeto
El objeto nulo o vacío, NULL

La mayoría de los objetos pueden tener atributos, como pueden ser las dimensiones o las etiquetas. De esta manera se pueden crear objetos con estructuras más complejas:

Matrices (resp. arrays), son vectores ligados a un atributo dim incluyendo las 2 (resp. -) dimensiones del objeto. Podemos crear matrices mediante la función matrix.
Factores, vectores con un atributo levels que corresponde a las posibles categorías de una variable categórica y un atributo class="factor". Por defecto son nominales, pero también pueden tener categorías ordenadas, caso en el que el atributo class es el vector c("ordered"," factor")
Fechas/Horas: hay 2 representaciones principales (ejecutando ?DateTimeClasses puede verse más información)
- vectores de tipo numérico (o entero), representando el número de segundos (resp. días) desde 1970-01-01 UTC, con clase POSIXct (resp. Date)
- listas, cuyos objetos son los distintos elementos del momento representado (año, mes, día, hora, minuto, segundo, etc.), con clase POSIXlt
Data frames, el tipo básico de base de datos en R. Es una lista de vectores (las columnas de la base de datos) -que pueden ser matrices o factores- que tienen todos la misma longitud, con un atributo names que designa los nombres de las columnas o variables de la base de datos y un atributo row.names que designa las filas u observaciones. Podemos crearlos mediante la función data.frame

Disponemos de un conjunto de operadores en R que nos permiten acceder y manipular los objetos. Además de los operadores aritméticos y lógicos, entre otros disponemos de los siguientes:

<-, para asignar un nombre a un objeto. Por ejemplo, x <- c(3,log(4), 5.2), y <- list(3, 'a', list(TRUE)), M <- matrix(c(1,2,3,4), nrow = 2, ncol = 2, byrow=TRUE)
[, para acceder a un elemento de un vector o a una sublista. Por ejemplo, x[3], y[2], y[3], M[,2] o M[1,2]
[[, para acceder a un elemento en una lista. Por ejemplo, y[[2]]

Otras características y algunas funciones prácticas

Los missings en R se representan con el símbolo NA, existente para los distintos tipos básicos, o también NaN en el caso de objetos numéricos.
Podemos nombrar los elementos de un vector o lista con la forma nombre = elemento. Por ejemplo y <- list(numero = 3, `letra` = 'a', "sublista" = list(TRUE)).
Cuando tenemos un vector o lista con nombres (como es el caso de cualquier data frame):
- podemos usar el nombre (como cadena de caracteres) en lugar del índice para acceder a los elementos del vector o a las sublistas. Por ejemplo y[c("letra", 'sublista')].
- en el caso de listas también podemos acceder a sus elementos mediante los operadores[[ o $. Por ejemplo y[['letra']] o y[["sublista"]], o equivalentemente y$letra o y$sublista (o y$`letra` o y$`sublista`).
Habitualmente es suficiente ejecutar el nombre de un objeto para ver el objeto. En ciertas ocasiones conviene evaluarlo explícitamente con la función print.
Otras funciones (principalmente informativas):
- Podemos informarnos acerca de la longitud, clase, estructura y atributos de un objeto por medio de las funciones length, class, str y attributes respectivamente.
- La función names nos permite conocer los nombres de los elementos de un objeto, como pueden ser las variables de un data frame.
- La función levels evaluada en un factor retorna las categorías del factor.
- Mediante dim vemos las dimensiones de un objeto, como pueden ser las filas y columnas de una matriz o de un data frame. Para objetos similares, dimnames informa acerca de los nombres de las filas y columnas.
Funciones de tipo is y as.
- Las funciones is permiten verificar que un objeto es de una determinada clase. Algunos ejemplos:
  - Es missing un elemento?: is.na(y)
  - Es un vector de tipo carácter o numérico?: is.character(y$numero), is.numeric(y$letra)
  - Es un data frame?: is.data.frame(y)
- Las funciones as fuerzan a los objetos a transformarse a la clase especificada cuando esto es posible.
  - Transformar a cadenas de caracteres. Por ejemplo, as.character(y)
  - Transformar a números. Por ejemplo, as.numeric(list("3", TRUE, FALSE, "a"))
  - Interpretar una cadena de texto como fecha. Por ejemplo, as.Date("01/05/1965", "%m/%d/%Y")
  - Transformar una matriz a data frame, por ejemplo as.data.frame(M)
Funciones para generar sucesiones de números y generar nuevos vectores:
- seq genera una sucesión de números desde uno inicial hasta uno final, y permite especificar o el tamaño de los intervalos o la longitud de la sucesión. Por ejemplo seq(from = 5, to = 20, by = 2)
- el operador : también permite generar sucesiones de números enteros. Por ejemplo, 5:20
- rep genera un vector de elementos a partir de un vector inicial dado. Por ejemplo, rep(x = c(1,3,5), times = 2) o rep(c(1,3,5), each = 5).
- paste y paste0 combinan cadenas de caracteres según se especifique. Por ejemplo, paste(y$letra, c(1,2))

Primeros pasos con una base de datos.

Asignamos la base de datos a una variable que llamamos dataw1. A partir de este momento tenemos un objeto de R con nombre dataw1, ya no es una base de datos en excel.

dataw1 <- read_excel("data_curs_stat/EP1.xls", sheet = 1) # aunque en este caso no es necesario, especificamos también que abrimos la hoja 1 del fichero excel
class(dataw1)

## [1] "tbl_df"     "tbl"        "data.frame"

length(dataw1)

## [1] 12

names(dataw1)

##  [1] "q0002_hhid"  "sex"         "age"         "mar_1"       "edu1"        "phys_hea1"   "hea1"        "dep1"        "score_lon1"  "score_sup1"  "income1"     "income_inf1"

str(dataw1)

## tibble [4,753 × 12] (S3: tbl_df/tbl/data.frame)
##  $ q0002_hhid : num [1:4753] 1 2 3 4 5 6 7 8 9 10 ...
##  $ sex        : chr [1:4753] "fem" "masc" "masc" "masc" ...
##  $ age        : chr [1:4753] "35-49" "35-49" "65-79" "50-64" ...
##  $ mar_1      : chr [1:4753] "Yes" "Yes" "Yes" "Yes" ...
##  $ edu1       : chr [1:4753] "Tertiary" "Tertiary" "Primary" "Primary" ...
##  $ phys_hea1  : chr [1:4753] "None" "None" "None" "1 physical health problem" ...
##  $ hea1       : num [1:4753] 70.5 78.9 66.6 77.4 52.5 ...
##  $ dep1       : chr [1:4753] "No" "No" "No" "No" ...
##  $ score_lon1 : num [1:4753] 3 3 3 3 3 3 6 3 3 6 ...
##  $ score_sup1 : num [1:4753] 12 12 12 11 12 13 10 12 11 13 ...
##  $ income1    : num [1:4753] 15 17 15 18 29 25 12 14 26 10 ...
##  $ income_inf1: chr [1:4753] "Normal" "Normal" "Normal" "Normal" ...

head(dataw1) # primeras (6) filas

## # A tibble: 6 × 12
##   q0002_hhid sex   age   mar_1 edu1     phys_hea1                  hea1 dep1  score_lon1 score_sup1 income1 income_inf1
##        <dbl> <chr> <chr> <chr> <chr>    <chr>                     <dbl> <chr>      <dbl>      <dbl>   <dbl> <chr>      
## 1          1 fem   35-49 Yes   Tertiary None                       70.5 No             3         12      15 Normal     
## 2          2 masc  35-49 Yes   Tertiary None                       78.9 No             3         12      17 Normal     
## 3          3 masc  65-79 Yes   Primary  None                       66.6 No             3         12      15 Normal     
## 4          4 masc  50-64 Yes   Primary  1 physical health problem  77.4 No             3         11      18 Normal     
## 5          5 fem   50-64 No    Tertiary None                       52.5 No             3         12      29 Normal     
## 6          6 fem   50-64 No    Tertiary None                       53.7 No             3         13      25 Normal

# head(dataw1, 3) # primeras 3 filas
# tail(dataw1) # últimas (6) filas
# View(dataw1) # vemos la base de datos completa

Acceso a las variables de una base de datos: los operadores `$` y `[[`

# dataw1$phys_hea1
# dataw1[["phys_hea1"]]
class(dataw1$phys_hea1)

## [1] "character"

#summary(dataw1$phys_hea1)
str(dataw1$phys_hea1) # chr quiere decir que phys_hea1 es una variable carácter, es decir, un vector de tipo cadena de carácteres

##  chr [1:4753] "None" "None" "None" "1 physical health problem" "None" "None" "None" "None" "None" "2+ physical health problems" "None" "None" "2+ physical health problems" "None" NA NA "2+ physical health problems" "None" ...

attributes(dataw1$phys_hea1)

## NULL

table(dataw1$phys_hea1)

## 
##   1 physical health problem 2+ physical health problems                        None 
##                        1538                         887                        1723

table(dataw1$phys_hea1, useNA = "ifany")

## 
##   1 physical health problem 2+ physical health problems                        None                        <NA> 
##                        1538                         887                        1723                         605

prop.table(table(dataw1$phys_hea1))

## 
##   1 physical health problem 2+ physical health problems                        None 
##                   0.3707811                   0.2138380                   0.4153809

Es posible introducir diferentes instrucciones en una misma línea separadas por ; (indiferentemente de los espacios en blanco en medio)

class(dataw1$hea1); str(dataw1$hea1) # num quiere decir que hea1 es una variable numérica

## [1] "numeric"

##  num [1:4753] 70.5 78.9 66.6 77.4 52.5 ...

summary(dataw1$hea1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   42.83   55.64   53.97   66.06   92.82     187

#?table

Otras funciones y operadores lógicos y aritméticos en R

Operadores aritméticos: +, -, *, /, %% (residuo), ^ (potencia)
Algunas funciones aritméticas: sum, round, floor, log, log10, log2, exp, sqrt, factorial
Operadores lógicos: == (igualdad), ! (negación), != (no igualdad), <, <=, >, >= (desigualdades), %in% (pertenencia), |, || (o), &, && (y)

3 %in% c(1:5)

## [1] TRUE

!3 %in% c(1:5)

## [1] FALSE

3 %in% c(5:50)

## [1] FALSE

!3 %in% c(5:50)

## [1] TRUE

TRUE & TRUE

## [1] TRUE

TRUE & FALSE

## [1] FALSE

FALSE | FALSE

## [1] FALSE

FALSE | TRUE

## [1] TRUE

x <- c(1,3,5,7, NA, 8,5, NA)
sum(is.na(x))

## [1] 2

Ejercicios

Instalar la librería openxlsx y buscar la ayuda de sus funciones read.xlsx y loadWorkbook
Buscar en la ayuda de las funciones presentadas (por ejemplo paste, seq o rep) qué opciones presentan según como se especifiquen sus argumentos. Para qué sirven las funciones seq_len y seq_along?
Usando las funciones setwd y read_excel convenientemente (hay varias posibilidades), guardar también las bases de datos EP2.xls y EP3.xls en dos variables que llamaremos dataw2 y dataw3. Podéis ver dataw1, dataw2, dataw3 y sus dimensiones en el panel Environment de RStudio?
Si hemos afirmado que un vector creado con la función c tiene todos los elementos del mismo tipo y que NA es un elemento de tipo lógico, por qué la instrucción c(1,3,5,7, NA, 8,5, NA) no da lugar a un error?
Los operadores lógicos | y || , y & y && se diferencian cuando comparan vectores lógicos con más de un elemento, como por ejemplo c(TRUE, FALSE). Cuáles son esas diferencias?

Nota: para buscar la ayuda de un operador como | podemos ejecutar la instrucción ?`|`

Transformación de variables. Factores

La función as.factor permite transformar a un factor, objeto que codifica una variable categórica:

class(as.factor(dataw1$phys_hea1))

## [1] "factor"

str(as.factor(dataw1$phys_hea1))

##  Factor w/ 3 levels "1 physical health problem",..: 3 3 3 1 3 3 3 3 3 2 ...

attributes(as.factor(dataw1$phys_hea1))

## $levels
## [1] "1 physical health problem"   "2+ physical health problems" "None"                       
## 
## $class
## [1] "factor"

table(as.factor(dataw1$phys_hea1), useNA = "ifany")

## 
##   1 physical health problem 2+ physical health problems                        None                        <NA> 
##                        1538                         887                        1723                         605

levels(as.factor(dataw1$phys_hea1))

## [1] "1 physical health problem"   "2+ physical health problems" "None"

Si nos interesa un atributo concreto, usamos la función attr:

attr(as.factor(dataw1$phys_hea1), "class")

## [1] "factor"

attr(as.factor(dataw1$phys_hea1), "levels")

## [1] "1 physical health problem"   "2+ physical health problems" "None"

Transformación de variables. Factores ordenados

Del mismo modo que podemos encontrarnos con variables categóricas ordenadas, como es el caso de las escalas de tipo Likert, los factores pueden ser ordenados. En este caso podemos utilizar la función as.ordered, aunque ésta puede ordenar las categorías erróneamente.

attributes(as.ordered(dataw1$phys_hea1))

## $levels
## [1] "1 physical health problem"   "2+ physical health problems" "None"                       
## 
## $class
## [1] "ordered" "factor"

str(as.ordered(dataw1$phys_hea1))

##  Ord.factor w/ 3 levels "1 physical health problem"<..: 3 3 3 1 3 3 3 3 3 2 ...

Ante ello, es preferible construír el factor especificando el orden que nos interesa:

x <- factor(dataw1$phys_hea1, ordered = TRUE, levels = c("None", "1 physical health problem", "2+ physical health problems"))
str(x)

##  Ord.factor w/ 3 levels "None"<"1 physical health problem"<..: 1 1 1 2 1 1 1 1 1 3 ...

attributes(x)

## $levels
## [1] "None"                        "1 physical health problem"   "2+ physical health problems"
## 
## $class
## [1] "ordered" "factor"

levels(x)

## [1] "None"                        "1 physical health problem"   "2+ physical health problems"

table(x, useNA = "ifany")

## x
##                        None   1 physical health problem 2+ physical health problems                        <NA> 
##                        1723                        1538                         887                         605

Transformación de una variable en una base de datos (con la librería `dplyr`)

Guardamos la transformación en la base de datos. Para ello usaremos la función mutate de la librería dplyr, que permite realizar varias transformaciones separadas por comas

library(dplyr)
#?mutate

dataw1 <- mutate(dataw1, phys_hea1 = as.factor(phys_hea1))

# es lo mismo que:

dataw1 <- dataw1 %>% mutate(phys_hea1 = as.factor(phys_hea1))

# podemos escrbir el argumento a la derecha de %>% en las líneas inferiores

dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
  mutate(phys_hea1 = as.factor(phys_hea1)) # transformamos la variable y la guardamos con el mismo nombre

# finalmente la base de datos queda guardada con el mismo nombre mediante la asignación inicial (dataw1 <- ...)

El primer argumento o valor de la función mutate es la base de datos o data frame.
Los siguientes argumentos son las transformaciones variable_output = transformación separadas por comas.
Usando el operador %>%, la función situada en el lado derecho (en este caso mutate) ya reconoce que su primer argumento es el valor situado en el lado izquierdo de %>% (en este caso dataw1)
Usamos la librería dplyr porque incluye otras funciones que facilitan realizar diversas operaciones y que veremos más adelante, de manera que usemos un lenguaje consistente. Sin embargo, podríamos escribir la instrucción especificada encima sin acudir a ninguna librería:

# Forma simple
dataw1$phys_hea1 <- as.factor(dataw1$phys_hea1)

# Usando la función transform
dataw1 <- transform(dataw1, phys_hea1 = as.factor(phys_hea1))

Y a partir de la versión 4.1.0, R dispone del nuevo operador |>:

dataw1 <- dataw1 |> transform(phys_hea1 = as.factor(phys_hea1))

Pipe

head(names(dataw1))

## [1] "q0002_hhid" "sex"        "age"        "mar_1"      "edu1"       "phys_hea1"

dataw1 %>% # la pipa más habitual: envía los datos de la izquierda a la función que sigue (a la derecha o abajo) y retorna su evaluación
  names() %>% # esta pipa está incluida en la librería dplyr, pero es original de la librería magrittr
  head()

## [1] "q0002_hhid" "sex"        "age"        "mar_1"      "edu1"       "phys_hea1"

library(magrittr)
# dataw1 %T>% # como la pipa anterior, pero retorna de nuevo los datos de la izquierda
#   View() %>% # sirve para usar los datos en más de una función, cuando las intermedias no retornan nada, como View, plot, ...
#   names() %>%
#   head()

dataw1 %$% # para funciones sin un argumento para data.frame's
  table(edu1, phys_hea1)

##               phys_hea1
## edu1           1 physical health problem 2+ physical health problems None
##   Less primary                       445                         396  304
##   Primary                            493                         276  471
##   Secondary                          406                         169  625
##   Tertiary                           192                          46  320

dataw1 %<>%
  mutate(phys_hea1 = as.factor(phys_hea1)) # la asignación está incluida en la pipa; 

# es lo mismo que
# dataw1 <- dataw1 %>% mutate(phys_hea1 = as.factor(phys_hea1))

Nota: la pipa %>% puede ser problemática con algunas funciones, como por ejemplo save. La pipa de R base |>, aunque no se puede utilizar para determinadas estructuras más complejas, es más segura. Por ejemplo, dataw1 |> save(file = "datos.Rdata") es exactamente lo mismo que save(dataw1, file = "datos.Rdata), pero dataw1 %>% save(file = "datos.Rdata") implica una instrucción ligeramente distinta.

Transformación de varias variables en una base de datos con la librería `dplyr`

Como la variable q0002_hhid es un id podría interesarnos que fuera de clase cadena

#class(dataw1$q0002_hhid)
#?as.character

Varias transformaciones seguidas las podemos evaluar:

Aplicando varias veces mutate:

dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
  mutate(phys_hea1 = as.factor(phys_hea1))  %>% # transformamos la variable phys_hea1 y la guardamos con el mismo nombre, y entonces
  mutate(dep1 = as.factor(dep1)) %>% # transformamos la variable dep1 y la guardamos con el mismo nombre, y entonces
  mutate(q0002_hhid = as.character(q0002_hhid)) # transformamos la variable q0002_hhid y la guardamos con el mismo nombre

Pero en lugar de eso, habitualmente se separan las instrucciones por comas en un sólo mutate:

dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
  mutate(phys_hea1 = as.factor(phys_hea1), dep1 = as.factor(dep1), q0002_hhid = as.character(q0002_hhid)) # transformamos las variables phys_hea1,  dep1 y q0002_hhid y las guardamos con los mismos nombres

O bien, cuando estamos aplicando la misma transformación (en este caso as.factor) a varias variables (y guardándolas con el mismo nombre), podemos usar across(vector de variables, función) dentro de mutate:

dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
  mutate(across(c(phys_hea1, dep1), as.factor), q0002_hhid = as.character(q0002_hhid)) # evaluamos as.factor a través de las variables phys_hea1 y dep1 y transformamos también number_id

Nota: Dentro de la función across, podemos seleccionar las variables que nos interesan de formas alternativas a escribir un vector de variables como c(phys_hea1, dep1). Veremos una introducción a estas alternativas en la siguiente diapositiva y en los ejercicios

Otras operaciones con funciones de la librería `dplyr`: selección de variables

Selección de variables

dataw1 %>% #partimos de la base de datos dataw1 y entonces
  select(q0002_hhid, hea1, dep1) %>% # mantenemos sólo las columnas q0002_hhid, hea1 y dep1, y entonces
  names() # miramos qué nombres tienen

## [1] "q0002_hhid" "hea1"       "dep1"

dataw1 %>%
  select(ends_with("1")) %>% # seleccionamos todas las variables cuyo nombre acaba en 1, y entonces
  names() # miramos qué nombres tienen

## [1] "mar_1"       "edu1"        "phys_hea1"   "hea1"        "dep1"        "score_lon1"  "score_sup1"  "income1"     "income_inf1"

Nota: Podemos seleccionar variables por índice y por nombre. Por ejemplo, para las 3 primeras variables escribiríamos select(1:3).

Qué ocurre si escribimos un - delante de end_with?
De qué otras funciones similares a ends_with disponemos?
Seleccionar todas las variables numéricas sin nombrarlas

Pull. Una función en dplyr para el operador $

dataw1 %>% pull(phys_hea1) %>% str()

##  Factor w/ 3 levels "1 physical health problem",..: 3 3 3 1 3 3 3 3 3 2 ...

Posibles combinaciones distintas de edad, educación y depresión

dataw1 %>% distinct(age, edu1, dep1)

## # A tibble: 59 × 3
##    age   edu1      dep1 
##    <chr> <chr>     <chr>
##  1 35-49 Tertiary  No   
##  2 65-79 Primary   No   
##  3 50-64 Primary   No   
##  4 50-64 Tertiary  No   
##  5 50-64 Secondary No   
##  6 80+   Primary   No   
##  7 18-34 Tertiary  No   
##  8 18-34 Secondary No   
##  9 65-79 Secondary <NA> 
## 10 50-64 Tertiary  <NA> 
## # … with 49 more rows

Otras operaciones con funciones de la librería `dplyr` y funciones de R para el cálculo de estadísticos

Ejemplo 1

dataw1 %>% #partimos de la base de datos dataw1 y entonces
  select(q0002_hhid, hea1, dep1) %>% # mantenemos sólo las columnas q0002_hhid, hea1 y dep1, y entonces
  filter(dep1 == "Yes") %>% # nos quedamos con las filas de quienes padecen depresión, y entonces
  arrange(desc(hea1)) %>% # ordenamos las filas por de mayor a menor valor de estado de salud, y entonces
  head(3) # nos quedamos con los 3 pacientes con depresión con mayor valor de Health state

## # A tibble: 3 × 3
##   q0002_hhid  hea1 dep1 
##        <dbl> <dbl> <chr>
## 1       4941  82.1 Yes  
## 2       4339  80.2 Yes  
## 3       1299  79.8 Yes

Ejemplo 2 (uso de funciones de R como mean, median o quantile; se pueden usar otras como var, sd, min, max, IQR, etc.)

dataw1 %>% #partimos de la base de datos dataw1 y entonces
  select(q0002_hhid, hea1, dep1, score_lon1) %>% # mantenemos sólo las columnas number_id, hea1 y dep1, y entonces
  group_by(dep1) %>% # agrupamos por las categorías de depresión, y entonces
  summarise(mean_hea1 = mean(hea1, na.rm = TRUE), median_hea1 = median(hea1, na.rm = TRUE), tercil_hea1 = quantile(hea1, probs = 1/3, na.rm = TRUE), n = n(), distinct_score_lon1 = n_distinct(score_lon1))

## # A tibble: 3 × 6
##   dep1  mean_hea1 median_hea1 tercil_hea1     n distinct_score_lon1
##   <chr>     <dbl>       <dbl>       <dbl> <int>               <int>
## 1 No         55.9        57.4        50.2  4073                   8
## 2 Yes        38.4        37.9        31.1   510                   8
## 3 <NA>      NaN          NA          NA     170                   1

# para cada categoría de depresión calculamos media, mediana y primer tercil de hea1,
# número de observaciones y número de observaciones distintas de score_lon1

dataw1 %>% distinct(dep1, score_lon1) %>% group_by(dep1) %>% summarise(n = n())

## # A tibble: 3 × 2
##   dep1      n
##   <chr> <int>
## 1 No        8
## 2 Yes       8
## 3 <NA>      1

Nota: usamos na.rm = TRUE dentro de mean, de median y de quantile para que calcule la media de aquellos valores que no son missing. En caso contrario, cuando hay missings el resultado es NA.

Nota: Todo lo anterior se puede realizar también con funciones de R sin necesidad de acudir a la librería dplyr, pero las alternativas pueden ser más complejas.

Cálculo de prevalencias:

Prevalencia de depresión:

dataw1 %>% 
  count(dep1) %>% # Contamos los individuos en cada categoría de depresión y entonces
  mutate(prop = 100*n/sum(n)) # calculamos el %

## # A tibble: 3 × 3
##   dep1      n  prop
##   <chr> <int> <dbl>
## 1 No     4073 85.7 
## 2 Yes     510 10.7 
## 3 <NA>    170  3.58

Prevalencia de depresión por grupos de edad y sin tener en cuenta los missings:

# ?is.na
dataw1 %>% 
  group_by(age) %>% 
  count(dep1) %>% # Contamos los individuos en cada categoría de depresión por grupo de edad y entonces
  filter(!is.na(dep1)) %>% # eliminamos los missings si nos interesa contar el porcentaje sobre el total de respuestas válidas, y entonces
  mutate(prop = 100*n/sum(n)) # Calculamos el %

## # A tibble: 10 × 4
## # Groups:   age [5]
##    age   dep1      n  prop
##    <chr> <chr> <int> <dbl>
##  1 18-34 No      382 93.9 
##  2 18-34 Yes      25  6.14
##  3 35-49 No      500 90.7 
##  4 35-49 Yes      51  9.26
##  5 50-64 No     1558 88.5 
##  6 50-64 Yes     202 11.5 
##  7 65-79 No     1298 87.3 
##  8 65-79 Yes     188 12.7 
##  9 80+   No      335 88.4 
## 10 80+   Yes      44 11.6

Vemos las variables y guardamos los resultados

Hemos dicho que R es como una calculadora, y que si no asignamos los objetos a variables, se muestran en consola pero no quedan guardados en ningún sitio.

Podemos ver los objetos que hemos creado en el panel Environment de RStudio o podemos ver sus nombres todavía con la función ls. Se pueden eliminar objetos mediante la función rm, separando sus nombres por comas.
Los objetos de R (como una base de datos o dataframe) se pueden guardar con la función save en ficheros con la extensión .rda o .rdata (aunque a veces también se escribe la r e incluso la d en mayúscula, por ejemplo .RData).
Los ficheros donde se guardan objetos de R se abren o cargan con la función load (o de forma distinta, con la función attach).
Los resultados, si son matrices, data frames o tablas, además, también se pueden guardar en ficheros .csv o con formatos de excel. De hecho, a menudo las funciones mencionadas de la forma read_ o read. tienen funciones correspondientes write_ o write.. Por ejemplo, podríamos guardar el data frame dataw1 en un fichero csv con write.csv:

# write.csv(dataw1, file = "datos_ola1.csv")

Los resultados y más generalmente todo aquello que aparece por consola puede ser guardados en ficheros de texto con la función sink. El primer argumento de esta función es el nombre del fichero donde se guardarán los resultados.
A no ser que se especifique la dirección del fichero, serán guardados en el directorio donde R está trabajando (getwd())

#?sink
sink("Summary_Estudi_poblacional_w1.txt", split = TRUE) # El fichero de texto donde se guardarán los resultados se llamará Summary_Estudi_poblacional_w1.txt, 
# pero podría llamarse como cualquier otro fichero de texto, no hay más restricción para el nombre que las del sistema operativo donde se está trabajando
table(dataw1$phys_hea1)

## 
##   1 physical health problem 2+ physical health problems                        None 
##                        1538                         887                        1723

summary(dataw1$hea1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   42.83   55.64   53.97   66.06   92.82     187

sink()

Ejercicios

Buscar todas las variables categóricas en las 3 bases de datos y transformar estas bases de datos de manera que esas variables sean factores. Mostrar las frecuencias de las categorías de algunas de esas variables. Transformar la variable q0002_hhid a cadena de carácteres.
Qué hace la función transmute?
Qué ocurre si escribimos un - delante de end_with en dataw1 %>% select(ends_with("1")) %>% names()?
De qué otras funciones similares a ends_with disponemos?
Seleccionar todas las variables numéricas sin nombrarlas
Comprender los distintos pasos en la instrucción dataw1 %>% distinct(dep1, score_lon1) %>% group_by(dep1) %>% summarise(n = n()) y las distintas funciones mencionadas referentes a estadísticos. Mirar también qué hace y cómo se puede usar la función range
Qué aparece en el fichero Summary_Estudi_poblacional_w1.txt (o como lo hayáis llamado)?

Condiciones y loops

Condiciones: if, else, ifelse (e if_else, la función equivalente en dplyr), case_when (otra función en dplyr)
Loops: for, while

y <- x <- as.character(dataw1$edu1); head(x, 10); table(x, useNA = "always")

##  [1] "Tertiary"  "Tertiary"  "Primary"   "Primary"   "Tertiary"  "Tertiary"  "Secondary" "Primary"   "Secondary" "Primary"

## x
## Less primary      Primary    Secondary     Tertiary         <NA> 
##         1364         1426         1329          629            5

for(i in seq_along(x)){
  if(is.na(x[i])){
    y[i] <- "Unknown"
  } else if (x[i] == "Less primary" || x[i] == "Primary"){
    y[i] <- "Until primary"
  } else if (x[i] == "Secondary" || x[i] == "Tertiary"){
    y[i] <- "Secondary or higher"
  } else {
    print(x[i]); y[i] <- x[i]
  }
}
head(y, 10); table(y, useNA = "always")

##  [1] "Secondary or higher" "Secondary or higher" "Until primary"       "Until primary"       "Secondary or higher" "Secondary or higher" "Secondary or higher" "Until primary"       "Secondary or higher" "Until primary"

## y
## Secondary or higher             Unknown       Until primary                <NA> 
##                1958                   5                2790                   0

z <- ifelse(x == "Secondary" | x == "Tertiary", "Secondary or higher", "Until primary")
head(z, 10); table(z, useNA = "always")

##  [1] "Secondary or higher" "Secondary or higher" "Until primary"       "Until primary"       "Secondary or higher" "Secondary or higher" "Secondary or higher" "Until primary"       "Secondary or higher" "Until primary"

## z
## Secondary or higher       Until primary                <NA> 
##                1958                2790                   5

w <- case_when(is.na(x) ~ "Unknown", 
               x %in% c("Less primary", "Primary") ~ "Until primary", 
               x %in% c("Secondary", "Tertiary") ~ "Secondary or higher", 
               TRUE ~ x)
head(w, 10); table(w, useNA = "always")

##  [1] "Secondary or higher" "Secondary or higher" "Until primary"       "Until primary"       "Secondary or higher" "Secondary or higher" "Secondary or higher" "Until primary"       "Secondary or higher" "Until primary"

## w
## Secondary or higher             Unknown       Until primary                <NA> 
##                1958                   5                2790                   0

Funciones *apply

Para aplicar una función a una matriz/data.frame/lista R dispone de las funciones de la familia apply:
- lapply: Devuelve una lista de resultados al aplicar una función a un vector o lista.
- sapply: Si es posible*, devuelve un vector de resultados al aplicar una función a un vector o lista.
- apply: Devuelve el vector o la lista de valores obtenidos aplicando una función a un vector o matriz por grupos.
- tapply: Devuelve una lista con los resultados obtenidos al aplicar una función a un vector por grupos.
- mapply: Devuelve un vector resultado de aplicar una función a tuplas de elementos en paralelo.

str(dataw1[,c("score_lon1", "score_sup1")])

## tibble [4,753 × 2] (S3: tbl_df/tbl/data.frame)
##  $ score_lon1: num [1:4753] 3 3 3 3 3 3 6 3 3 6 ...
##  $ score_sup1: num [1:4753] 12 12 12 11 12 13 10 12 11 13 ...

lapply(dataw1[,c("score_lon1", "score_sup1")], class)

## $score_lon1
## [1] "numeric"
## 
## $score_sup1
## [1] "numeric"

sapply(dataw1[,c("score_lon1", "score_sup1")], class)

## score_lon1 score_sup1 
##  "numeric"  "numeric"

str(apply(dataw1[, c("score_lon1", "score_sup1")], 1, sum))

##  num [1:4753] 15 15 15 14 15 16 16 15 14 19 ...

apply(dataw1[, c("score_lon1", "score_sup1")], 2, sum, na.rm = TRUE)

## score_lon1 score_sup1 
##      16712      50008

mapply(paste, 1:4, 4:1)

## [1] "1 4" "2 3" "3 2" "4 1"

* Del mismo modo que la función c, sapply intenta forzar a todos los elementos de la lista inicialmente resultante a tener el mismo tipo, en el orden NULL < lógico < entero < numérico < carácter (esencialmente)

Variables y dimensiones de las bbdd

names(dataw1)

##  [1] "q0002_hhid"  "sex"         "age"         "mar_1"       "edu1"        "phys_hea1"   "hea1"        "dep1"        "score_lon1"  "score_sup1"  "income1"     "income_inf1"

dim(dataw1)

## [1] 4753   12

ncol(dataw1)

## [1] 12

nrow(dataw1)

## [1] 4753

names(dataw2)

## [1] "q0002_hhid" "dep2"       "score_lon2" "score_sup2" "income2"    "phys_hea2"  "hea2"

dim(dataw2)

## [1] 2400    7

names(dataw3)

## [1] "q0002_hhid" "dep3"       "score_lon3" "score_sup3" "income3"    "phys_hea3"  "health3"

dim(dataw3)

## [1] 1512    7

Fusión (merge) de bbdd

Cuántos missings tiene la variable q0002_hhid?

table(is.na(dataw1$q0002_hhid)); table(is.na(dataw2$q0002_hhid)); table(is.na(dataw3$q0002_hhid))

## 
## FALSE 
##  4753

## 
## FALSE 
##  2400

## 
## FALSE 
##  1512

Ya tenemos las bbdd abiertas y guardadas. Sólo tenemos que unir por la(s) variable(s) identificadora(s), en este caso q0002_hhid.

data <- dataw1 %>% # partimos de la base de datos dataw1 y entonces
  full_join(dataw2, by = "q0002_hhid") %>% # unimos horizontalmente con todas las observaciones de dataw2 con q0002_hhid iguales a los de dataw1 y añadimos las nuevas, y entonces
  full_join(dataw3, by = "q0002_hhid") # unimos horizontalmente con todas las observaciones de dataw3 con q0002_hhid iguales a los que ya había y añadimos las nuevas
names(data)

##  [1] "q0002_hhid"  "sex"         "age"         "mar_1"       "edu1"        "phys_hea1"   "hea1"        "dep1"        "score_lon1"  "score_sup1"  "income1"     "income_inf1" "dep2"        "score_lon2"  "score_sup2"  "income2"    
## [17] "phys_hea2"   "hea2"        "dep3"        "score_lon3"  "score_sup3"  "income3"     "phys_hea3"   "health3"

dim(data)

## [1] 4753   24

?full_join

Bases de datos en formato de Stata

Abrimos las 3 olas de un ensayo clínico en fichero separados. Usamos la función read_dta en la librería haven.

Puede abrir bases de datos guardadas con diferentes versiones de Stata, pero si fueron guardadas con versiones previas a la 14 en Mac o Linux, entonces hay que ajustar el argumento encoding = "latin1" (información que proporciona con más detalle la ayuda de la función, ?read_dta)
La función write_dta permite guardar data frames en formato de Stata para distintas versiones, especificando el argumento version.
Otras librerías que permiten abrir ficheros de Stata son foreign y readstata13.

library(haven)
ac1 <- read_dta("data_curs_stat/Assaig_clinic_w1.dta"); ac2 <- read_dta("data_curs_stat/Assaig_clinic_w2.dta"); ac3 <- read_dta("data_curs_stat/Assaig_clinic_w3.dta")

Una base de datos en formato de Stata está guardada de manera distinta que una base de datos en formato de excel. Por ello, al abrirla como un data frame, éste tiene una estructura ligeramente distinta.

Recordemos que cuando abrimos la base de datos que guardamos en el objeto dataw1, la variable phys_hea1 era un vector carácter y sus elementos eran las distintas categorías, de manera que la transformamos de manera razonable en factor mediante la función as.factor. Ahora, sin embargo, al abrir bases de datos en formato Stata con la librería haven encontramos algunas diferencias

class(ac1$edu1)

## [1] "haven_labelled" "vctrs_vctr"     "double"

str(ac1$edu1)

##  dbl+lbl [1:400] 2, 1, 1, 3, 1, 2, 3, 2, 1, 1, 3, 2, 2, 1, 1, 3, 2, 1, 2, 2, 1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 3, 1, 3, 1, 1, 2, 4, 3, 4, 3, 1, 2, 1, 2, 3, 1, 3, 1, 2, 2, 1, 1, 3, 1, 1, 3, 1, 1, 1, 2, 1, 2, 2, 1, 1, 2, 3, 2, 1, 1, 1, 1, 1,...
##  @ label       : chr "Level of education"
##  @ format.stata: chr "%12.0g"
##  @ labels      : Named num [1:4] 1 2 3 4
##   ..- attr(*, "names")= chr [1:4] "Less primary" "Primary" "Secondary" "Tertiary"

attributes(ac1$edu1)

## $label
## [1] "Level of education"
## 
## $format.stata
## [1] "%12.0g"
## 
## $class
## [1] "haven_labelled" "vctrs_vctr"     "double"        
## 
## $labels
## Less primary      Primary    Secondary     Tertiary 
##            1            2            3            4

table(ac1$edu1)

## 
##   1   2   3   4 
## 179 112  89  19

Variables con etiquetas

En formatos de fichero distintos de excel o .csv, como pueden ser los de Stata o SPSS, las categorías de las variables categóricas pueden estar grabadas como etiquetas, mientras los valores son numéricos.

class(as.factor(ac1$edu1))

## [1] "factor"

str(as.factor(ac1$edu1))

##  Factor w/ 4 levels "1","2","3","4": 2 1 1 3 1 2 3 2 1 1 ...

attributes(as.factor(ac1$edu1))

## $levels
## [1] "1" "2" "3" "4"
## 
## $class
## [1] "factor"

table(as.factor(ac1$edu1))

## 
##   1   2   3   4 
## 179 112  89  19

En estos casos, la función as_factor de la librería haven es más adecuada que as.factor para recuperar las etiquetas, ya que, a diferencia de esta última, tiene en cuenta los atributos de la variable original.

class(as_factor(ac1$edu1))

## [1] "factor"

str(as_factor(ac1$edu1))

##  Factor w/ 4 levels "Less primary",..: 2 1 1 3 1 2 3 2 1 1 ...
##  - attr(*, "label")= chr "Level of education"

attributes(as_factor(ac1$edu1))

## $levels
## [1] "Less primary" "Primary"      "Secondary"    "Tertiary"    
## 
## $class
## [1] "factor"
## 
## $label
## [1] "Level of education"

table(as_factor(ac1$edu1))

## 
## Less primary      Primary    Secondary     Tertiary 
##          179          112           89           19

Combinar bases de datos

dim(ac1); names(ac1); dim(ac2); names(ac2); dim(ac3); names(ac3)

## [1] 400   9

## [1] "q0002_hhid" "number_id"  "ID_ECS"     "sex"        "age"        "mar1"       "edu1"       "hea1"       "grups"

## [1] 400   6

## [1] "q0002_hhid" "number_id"  "ID_ECS"     "grups"      "hea2"       "dep2"

## [1] 400   6

## [1] "q0002_hhid" "number_id"  "ID_ECS"     "grups"      "hea3"       "dep3"

Para poder combinarlas verticalmente, los nombres de las columnas que queremos combinar tienen que ser iguales. Para ello usamos otra función de la librería dplyr: rename

acr1 <- ac1 %>% # partimos de ac1 y entonces
  rename(hea = hea1) # renombramos hea1 como hea
acr2 <- ac2 %>% # partimos de ac2 y entonces
  rename(hea = hea2, dep = dep2) # renombramos hea2 como hea y dep2 como dea
acr3 <- ac3 %>% 
  rename(hea = hea3, dep = dep3) # mismo proceso para ac3

#names(ac1); names(ac2); names(ac3)

Finalmente combinamos las filas con otra función de la librería dplyr: bind_rows

#?bind_rows
acv <- bind_rows(acr1, acr2, acr3, .id = "wave")

Otras opciones para combinar bases de datos

Nota: La función rbind de R hace esencialmente lo mismo, pero necesita que las 3 bases de datos tengan exactamente las mismas columnas. En caso en que esto no ocurre, como el presente, da un error. Además, tanto en el caso de rbind como de bind_rows, conviene que las columnas por las que se combina (las de igual nombre) tengan la misma clase, pues en caso contrario pueden ocurrir errores o comportamiento extraños.

Nota: Las funciones cbind y bind_cols de R y dplyr respectivamente combinan por columnas. Se diferencian de un merge en que no hay una columna “común” por la que fusionar, sino que se añaden las columnas de los distintos objetos, tal como están ordenadas en cada uno de ellos. Además, todas las columnas tienen que tener el mismo número de elementos.

Fusionamos horizontalmente las 3 bbdd originales de los ensayos clínicos. En este caso no las podemos combinar porque tienen distinto número de filas:

ach <- ac1 %>%
  full_join(ac2, by = c("q0002_hhid", "number_id", "ID_ECS", "grups")) %>%
  full_join(ac3, by = c("q0002_hhid", "number_id", "ID_ECS", "grups"))
ach %>% 
  head()

## # A tibble: 6 × 13
##   q0002_hhid number_id ID_ECS       sex       age      mar1             edu1  hea1            grups  hea2      dep2  hea3      dep3
##        <dbl> <chr>      <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl>        <dbl+lbl> <dbl>        <dbl+lbl> <dbl> <dbl+lbl> <dbl> <dbl+lbl>
## 1       2306 724002306   2306   2 [fem] 3 [50-64]   0 [No]  2 [Primary]       45.9 1 [Intervention]  NA     NA       NA    NA      
## 2       3363 724003363   3363   2 [fem] 4 [65-79]   0 [No]  1 [Less primary]  23.3 1 [Intervention]  24.6    0 [No]  NA    NA      
## 3       3228 724003228   3228   2 [fem] 3 [50-64]   1 [Yes] 1 [Less primary]  49.1 1 [Intervention]  35.2    0 [No]  NA    NA      
## 4       2371 724002371   2371   2 [fem] 1 [18-34]   0 [No]  3 [Secondary]     57.5 1 [Intervention]  NA     NA       NA    NA      
## 5       2452 724002452   2452   2 [fem] 2 [35-49]   0 [No]  1 [Less primary]  41.3 1 [Intervention]  44.3    0 [No]  31.6   1 [Yes]
## 6       2750 724002750   2750   2 [fem] 5 [80+]     0 [No]  2 [Primary]       19.5 1 [Intervention]  NA     NA       NA    NA

cbind(ac1, ac2, ac3) %>% 
  head()

##   q0002_hhid number_id ID_ECS sex age mar1 edu1     hea1 grups q0002_hhid number_id ID_ECS grups     hea2 dep2 q0002_hhid number_id ID_ECS grups     hea3 dep3
## 1       2306 724002306   2306   2   3    0    2 45.88945     1         48 724000048     48     0 46.31060    0         48 724000048     48     0 47.35691    0
## 2       3363 724003363   3363   2   4    0    1 23.27672     1         61 724000061     61     0 53.24683    0         61 724000061     61     0 67.96230    0
## 3       3228 724003228   3228   2   3    1    1 49.08676     1         90 724000090     90     1       NA   NA         90 724000090     90     1       NA   NA
## 4       2371 724002371   2371   2   1    0    3 57.50416     1        128 724000128    128     0 60.88645    0        128 724000128    128     0 47.54218    0
## 5       2452 724002452   2452   2   2    0    1 41.27082     1        181 724000181    181     1       NA   NA        181 724000181    181     1       NA   NA
## 6       2750 724002750   2750   2   5    0    2 19.50435     1        191 724000191    191     0       NA   NA        191 724000191    191     0       NA   NA

Las funciones de las librerías de la familia tidyverse tratan los nombres de las variables de un data frame de una manera distinta a como lo hacen las funciones de R base. Por ejemplo, admiten determinados nombres especiales, pero sin embargo no admiten variables distintas con el mismo nombre, como sí ocurre en el último ejemplo.

Rellenando datos

Cuando se combinan las filas mediante bind_cols, las columnas que sólo están en una base de datos se llenan con missings. Por ejemplo, es el caso de las variables sociodemográficas, que en el dataframe acv sólo tienen datos en las filas correspondientes a la ola 1, que son las filas obtenidas de ac1.

acv %>%
  arrange(q0002_hhid) %>% # ordenamos por id
  head()

## # A tibble: 6 × 11
##   wave  q0002_hhid number_id ID_ECS       sex        age      mar1           edu1   hea       grups       dep
##   <chr>      <dbl> <chr>      <dbl> <dbl+lbl>  <dbl+lbl> <dbl+lbl>      <dbl+lbl> <dbl>   <dbl+lbl> <dbl+lbl>
## 1 1             48 724000048     48   2 [fem]  3 [50-64]    0 [No]  3 [Secondary]  47.9 0 [Control]   NA     
## 2 2             48 724000048     48  NA       NA           NA      NA              46.3 0 [Control]    0 [No]
## 3 3             48 724000048     48  NA       NA           NA      NA              47.4 0 [Control]    0 [No]
## 4 1             61 724000061     61   2 [fem]  2 [35-49]    0 [No]  4 [Tertiary]   63.2 0 [Control]   NA     
## 5 2             61 724000061     61  NA       NA           NA      NA              53.2 0 [Control]    0 [No]
## 6 3             61 724000061     61  NA       NA           NA      NA              68.0 0 [Control]    0 [No]

Acudimos ahora a la función fill de la librería tidyr:

library(tidyr)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:magrittr':
## 
##     extract

acv %>%
  group_by(q0002_hhid) %>% # agrupamos por id y entonces
  arrange(wave) %>% # ordenamos por ola para tener la fila con datos antes que las otras, y entonces
  fill(sex, age, mar1, edu1) %>% # para cada id rellenamos las filas vacías de las variables especificadas con los valores que no son missings, y entonces
  ungroup() %>% # desagrupamos y entonces
  arrange(q0002_hhid) %>% # ordenamos por id para ver las mismas filas que arriba, y entonces
  head() # nos quedamos con las primeras filas

## # A tibble: 6 × 11
##   wave  q0002_hhid number_id ID_ECS       sex       age      mar1          edu1   hea       grups       dep
##   <chr>      <dbl> <chr>      <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl>     <dbl+lbl> <dbl>   <dbl+lbl> <dbl+lbl>
## 1 1             48 724000048     48   2 [fem] 3 [50-64]    0 [No] 3 [Secondary]  47.9 0 [Control]   NA     
## 2 2             48 724000048     48   2 [fem] 3 [50-64]    0 [No] 3 [Secondary]  46.3 0 [Control]    0 [No]
## 3 3             48 724000048     48   2 [fem] 3 [50-64]    0 [No] 3 [Secondary]  47.4 0 [Control]    0 [No]
## 4 1             61 724000061     61   2 [fem] 2 [35-49]    0 [No] 4 [Tertiary]   63.2 0 [Control]   NA     
## 5 2             61 724000061     61   2 [fem] 2 [35-49]    0 [No] 4 [Tertiary]   53.2 0 [Control]    0 [No]
## 6 3             61 724000061     61   2 [fem] 2 [35-49]    0 [No] 4 [Tertiary]   68.0 0 [Control]    0 [No]

Ejercicios

Con full_join creamos una base de datos resultado de fusionar las 3 iniciales e incluír todas las observaciones de cada una de ellas. Consultando en la ayuda, este ejercicio consiste en fusionar las 3 bases de datos, pero incluyendo sólo aquellas observaciones comunes a las 3. Cuántas observaciones tiene?
Ver qué variables tiene acv; para qué se añadió .id = "wave"?; ejecutar rbind(ac1,ac2,ac3)
Sustituír cbind por bind_cols en el último ejemplo. Qué ocurre si se sustituye el valor por defecto del argumento .name_repair por los otros valores posibles?
Modificar acv de tal manera que sus variables sociodemográficas estén completas tal como se indicó en la diapositiva anterior.
Qué variables tiene acv?
Qué variables de acv no tienen necesariamente los mismos valores en las distintas olas?

Pivotar (vertical/largo –> horizontal/ancho)

Bases de datos horizontales/anchas: aquellas con una fila para cada individuo y tantas columnas como observaciones por variable, por ejemplo el producto de un una fusión, como es el caso de ach
Bases de datos verticales/largas (longitudinales): aquellas con tantas filas como observaciones se hayan hecho a todos los individuos y una columna por variable, por ejemplo el resultado de combinar filas, como es el caso de acv.

Pivotar es el proceso de transformar una base de datos vertical/longitudinal a una horizontal o a la inversa.

Pivotaje usando la función pivot_wider de la librería tidyr:

#?pivot_wider
acv %>%
  pivot_wider(names_from = c("wave"), values_from = c("hea", "dep")) %>%
  # select(-where(~all(is.na(.)))) %>% # quitamos aquellas columnas que cumplen que todos sus valores son missings
  head()

## # A tibble: 6 × 14
## # Groups:   q0002_hhid [6]
##   q0002_hhid number_id ID_ECS       sex       age      mar1             edu1            grups hea_1 hea_2 hea_3     dep_1     dep_2    dep_3
##        <dbl> <chr>      <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl>        <dbl+lbl>        <dbl+lbl> <dbl> <dbl> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lb>
## 1       2306 724002306   2306   2 [fem] 3 [50-64]   0 [No]  2 [Primary]      1 [Intervention]  45.9  NA    NA          NA   NA      NA      
## 2       3363 724003363   3363   2 [fem] 4 [65-79]   0 [No]  1 [Less primary] 1 [Intervention]  23.3  24.6  NA          NA    0 [No] NA      
## 3       3228 724003228   3228   2 [fem] 3 [50-64]   1 [Yes] 1 [Less primary] 1 [Intervention]  49.1  35.2  NA          NA    0 [No] NA      
## 4       2371 724002371   2371   2 [fem] 1 [18-34]   0 [No]  3 [Secondary]    1 [Intervention]  57.5  NA    NA          NA   NA      NA      
## 5       2452 724002452   2452   2 [fem] 2 [35-49]   0 [No]  1 [Less primary] 1 [Intervention]  41.3  44.3  31.6        NA    0 [No]  1 [Yes]
## 6       2750 724002750   2750   2 [fem] 5 [80+]     0 [No]  2 [Primary]      1 [Intervention]  19.5  NA    NA          NA   NA      NA

Pivotar (horizontal/ancho –> vertical/largo)

Si tenemos un sólo conjunto de columnas longitudinal, por ejemplo hea1, hea2 y hea3, podemos acudir a la función pivot_longer de la librería tidyr:

ach %>%
  select(-starts_with("dep")) %>% # quitamos las columnas que empiezan con dep y entonces
  pivot_longer(cols = starts_with("hea"), names_to = "wave", values_to = "hea") %>% # pivotamos, y entonces
  arrange(q0002_hhid) %>% head()

## # A tibble: 6 × 10
##   q0002_hhid number_id ID_ECS       sex       age      mar1          edu1       grups wave    hea
##        <dbl> <chr>      <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl>     <dbl+lbl>   <dbl+lbl> <chr> <dbl>
## 1         48 724000048     48   2 [fem] 3 [50-64]    0 [No] 3 [Secondary] 0 [Control] hea1   47.9
## 2         48 724000048     48   2 [fem] 3 [50-64]    0 [No] 3 [Secondary] 0 [Control] hea2   46.3
## 3         48 724000048     48   2 [fem] 3 [50-64]    0 [No] 3 [Secondary] 0 [Control] hea3   47.4
## 4         61 724000061     61   2 [fem] 2 [35-49]    0 [No] 4 [Tertiary]  0 [Control] hea1   63.2
## 5         61 724000061     61   2 [fem] 2 [35-49]    0 [No] 4 [Tertiary]  0 [Control] hea2   53.2
## 6         61 724000061     61   2 [fem] 2 [35-49]    0 [No] 4 [Tertiary]  0 [Control] hea3   68.0

Cuando hay más de un conjunto de columnas que queremos pivotar, como en este caso, en el que tenemos por un lado hea1, hea2 y hea3 y por otro lado, dep2 y dep3, usaremos la función reshape de R. Para ello, primero necesitamos el mismo número de columnas de hea y dep, por tanto creamos la columna dep1.

ach$dep1 <- NA  # como en un data frame todas las variables tienen las mismas filas, R ya entiende que tiene que asignar NA para cada fila.

También necesitamos que el objeto sea un data frame propiamente dicho, en lugar de un tibble (son así llamados los data frames creados por las funciones de la familia tidyverse, que se caracterizan porque su clase es c("tbl_df", "tbl", "data.frame"))

class(ach)

## [1] "tbl_df"     "tbl"        "data.frame"

class(as.data.frame(ach))

## [1] "data.frame"

reshape(as.data.frame(ach), varying = list(paste0("hea", 1:3), paste0("dep", 1:3)), v.names = c("hea", "dep"), timevar = "wave", idvar = "q0002_hhid", direction = "long") %>% 
  arrange(q0002_hhid) %>%
  head()

##      q0002_hhid number_id ID_ECS sex age mar1 edu1 grups wave      hea dep
## 48.1         48 724000048     48   2   3    0    3     0    1 47.94347  NA
## 48.2         48 724000048     48   2   3    0    3     0    2 46.31060   0
## 48.3         48 724000048     48   2   3    0    3     0    3 47.35691   0
## 61.1         61 724000061     61   2   2    0    4     0    1 63.24433  NA
## 61.2         61 724000061     61   2   2    0    4     0    2 53.24683   0
## 61.3         61 724000061     61   2   2    0    4     0    3 67.96230   0

Correlación y tests de normalidad

Cálculo de correlaciones

#cor(dataw1$score_sup, dataw1$income)
cor.test(dataw1$score_sup1, dataw1$income1)

## 
##  Pearson's product-moment correlation
## 
## data:  dataw1$score_sup1 and dataw1$income1
## t = 3.2915, df = 3518, p-value = 0.001007
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.02241180 0.08828366
## sample estimates:
##        cor 
## 0.05540802

Test de Shapiro-Wilk

shapiro.test(dataw1$hea1)

## 
##  Shapiro-Wilk normality test
## 
## data:  dataw1$hea1
## W = 0.98634, p-value < 2.2e-16

Test de Kolmogorov-Smirnov

ks.test(dataw1$hea1, "pnorm")

## Warning in ks.test(dataw1$hea1, "pnorm"): ties should not be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  dataw1$hea1
## D = 0.99954, p-value < 2.2e-16
## alternative hypothesis: two-sided

Ambos test son significativos, así pues descartamos la hipótesis nula consistente en la normalidad de la distribución de hea1 para la muestra en dataw1.
Para otros tests de normalidad como los de Lilliefors, Anderson-Darling o Jarque-Bera, disponemos de funciones en librerías como DescTools, nortest o tseries.

Tests chi-cuadrado y de Fisher

table(dataw1$mar_1, dataw1$sex)

##      
##        fem masc
##   No  1282  626
##   Yes 1320 1525

Test de chi (o ji o khi)-cuadrado de Pearson
- Prueba no paramétrica que se utiliza para probar la independencia de dos variables categóricas nominales entre sí mediante la presentación de los datos en tablas de contingencia.
- Hipótesis nula: independientes vs Hipotesis alternativa ($p\approx0$): dependientes
- En la tabla de las frecuencias esperadas cada celda debe ser mayor o igual que 5 para poder aplicar este test

test <- chisq.test(dataw1$mar_1, dataw1$sex)
test  # por defecto, test con correccion de continuidad si es 2x2

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  dataw1$mar_1 and dataw1$sex
## X-squared = 198.48, df = 1, p-value < 2.2e-16

test$expected # frecuencia esperada de cada celda

##             dataw1$sex
## dataw1$mar_1      fem      masc
##          No  1044.523  863.4774
##          Yes 1557.477 1287.5226

chisq.test(dataw1$mar_1, dataw1$sex, correct = FALSE) # test sin corrección

## 
##  Pearson's Chi-squared test
## 
## data:  dataw1$mar_1 and dataw1$sex
## X-squared = 199.31, df = 1, p-value < 2.2e-16

Test de Fisher

fisher.test(dataw1$mar_1, dataw1$sex)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  dataw1$mar_1 and dataw1$sex
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.093293 2.674704
## sample estimates:
## odds ratio 
##   2.365541

Tests univariantes

Test t de Student de comparación de una media muestral a la hipótesis nula

dataw1$new_var <- rnorm(nrow(dataw1))
mean(dataw1$new_var)

## [1] 0.00137449

t.test(dataw1$new_var, mu = 0)

## 
##  One Sample t-test
## 
## data:  dataw1$new_var
## t = 0.092092, df = 4752, p-value = 0.9266
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.02788569  0.03063467
## sample estimates:
##  mean of x 
## 0.00137449

t.test(dataw1$new_var)$conf.int # intervalo de confianza para la media muestral de new_var en dataw1

## [1] -0.02788569  0.03063467
## attr(,"conf.level")
## [1] 0.95

Test no paramétrico de los rangos con signo de Wilcoxon, de simetría de los datos respecto a la hipótesis nula

mean(dataw1$hea1, na.rm = TRUE)

## [1] 53.96553

wilcox.test(dataw1$hea1, mu = 54, conf.int = TRUE, conf.level = 0.95, exact = FALSE)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  dataw1$hea1
## V = 5409216, p-value = 0.0278
## alternative hypothesis: true location is not equal to 54
## 95 percent confidence interval:
##  54.06244 55.07854
## sample estimates:
## (pseudo)median 
##       54.57055

Como conclusiones, no podemos descartar que la media de new_var sea 0, pero sí podemos descartar que los valores de hea1 se localizen en torno a 54.

Prueba t de Student (o t-test)

Comparamos entre grupos mediante un t-test:

t.test(hea1 ~ sex, data = dataw1)

## 
##  Welch Two Sample t-test
## 
## data:  hea1 by sex
## t = -14.557, df = 4523.6, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group fem and group masc is not equal to 0
## 95 percent confidence interval:
##  -7.861186 -5.995119
## sample estimates:
##  mean in group fem mean in group masc 
##           50.82465           57.75280

t.test(data$hea1, data$hea2, paired = TRUE) # comparamos el estado de salud en la ola 1 con la ola 2

## 
##  Paired t-test
## 
## data:  data$hea1 and data$hea2
## t = 0.81332, df = 2384, p-value = 0.4161
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2407952  0.5820965
## sample estimates:
## mean of the differences 
##               0.1706506

t.test(hea1 ~ sex, data = dataw1, mu = -7) # comparamos la diferencia de medias con el valor mu (hipótesis nula)

## 
##  Welch Two Sample t-test
## 
## data:  hea1 by sex
## t = 0.15097, df = 4523.6, p-value = 0.88
## alternative hypothesis: true difference in means between group fem and group masc is not equal to -7
## 95 percent confidence interval:
##  -7.861186 -5.995119
## sample estimates:
##  mean in group fem mean in group masc 
##           50.82465           57.75280

Nota: mediante el argumento alternative de la función t.test, podemos especificar si la hipótesis nula es la igualdad (test bilateral) o una desigualdad.

Prueba Wilcoxon-Mann-Whitney (también test U de Mann-Whitney o test de la suma de rangos de Wilcoxon)

Prueba no paramétrica que compara los valores en dos grupos

wilcox.test(hea1 ~ sex, data = dataw1, 
            alternative = "two.sided", mu = 0, 
            paired = FALSE, exact = FALSE, 
            conf.int = TRUE, conf.level = 0.95)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  hea1 by sex
## W = 1964188, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -7.958215 -6.019921
## sample estimates:
## difference in location 
##              -6.985333

library(ggplot2)

Descartamos la hipótesis nula, por tanto podemos concluír que la distribución de hea1 en los hombres es significativamente diferente de la distribución en las mujeres.

dataw1 %>% 
  filter(!is.na(hea1)) %>% 
  ggplot(aes(x = hea1)) + # asignamos la variable hea1 al eje x
  geom_histogram(bins = 15, fill = "turquoise", color = "black") + # creamos el histograma y lo personalizamos 
  facet_wrap(~sex) # dividimos la gráfica para las distintas categorías de sex

dataw1 %>% 
  filter(!is.na(hea1)) %>% 
  ggplot(aes(x = sex, y = hea1)) +  # asignamos la variable sex al eje x y la variable hea1 al eje y
  geom_boxplot() + # creamos el boxplot
  coord_flip()

ANOVA

Test de comparación de la distribución de una variable normal entre múltiples grupos

ks.test(dataw1$new_var, "pnorm")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  dataw1$new_var
## D = 0.013135, p-value = 0.3851
## alternative hypothesis: two-sided

dataw1 %>%
  filter(!is.na(edu1)) %>%
  ggplot(aes(x = edu1, y = new_var)) +
  geom_boxplot()

anovafit <- aov(new_var ~ edu1, data = dataw1)
summary(anovafit)

##               Df Sum Sq Mean Sq F value Pr(>F)
## edu1           3      2   0.631   0.596  0.618
## Residuals   4744   5026   1.059               
## 5 observations deleted due to missingness

Regresión lineal

#?lm
fit <- lm(hea1 ~ age + sex + income1 + dep1, data = dataw1)
summary(fit)

## 
## Call:
## lm(formula = hea1 ~ age + sex + income1 + dep1, data = dataw1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.560  -8.844   0.710   9.444  40.641 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  64.99142    0.78354  82.946  < 2e-16 ***
## age35-49     -5.41946    0.95045  -5.702 1.28e-08 ***
## age50-64    -11.26679    0.79368 -14.196  < 2e-16 ***
## age65-79    -19.19973    0.80852 -23.747  < 2e-16 ***
## age80+      -30.15751    1.02706 -29.363  < 2e-16 ***
## sexmasc       5.07076    0.44117  11.494  < 2e-16 ***
## income1       0.16659    0.02742   6.075 1.37e-09 ***
## dep1Yes     -14.52514    0.67309 -21.580  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.09 on 3667 degrees of freedom
##   (1078 observations deleted due to missingness)
## Multiple R-squared:  0.3837, Adjusted R-squared:  0.3825 
## F-statistic: 326.1 on 7 and 3667 DF,  p-value: < 2.2e-16

coef(fit) # o fit$coefficients

## (Intercept)    age35-49    age50-64    age65-79      age80+     sexmasc     income1     dep1Yes 
##   64.991424   -5.419455  -11.266791  -19.199732  -30.157514    5.070758    0.166592  -14.525141

#?confint
confint(fit, level = 0.95) # intervalos de confianza

##                   2.5 %      97.5 %
## (Intercept)  63.4552039  66.5276443
## age35-49     -7.2829188  -3.5559915
## age50-64    -12.8228922  -9.7106894
## age65-79    -20.7849349 -17.6145292
## age80+      -32.1711810 -28.1438463
## sexmasc       4.2058013   5.9357156
## income1       0.1128245   0.2203594
## dep1Yes     -15.8448166 -13.2054651

#predict
plot(fit)

Regresión logística

#?glm
fit <- glm(dep1 ~ age + sex + income1*hea1 + mar_1*score_lon1, data = dataw1, family = binomial)
summary(fit)

## 
## Call:
## glm(formula = dep1 ~ age + sex + income1 * hea1 + mar_1 * score_lon1, 
##     family = binomial, data = dataw1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0231  -0.4485  -0.2799  -0.1678   3.2085  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.9027594  0.4760209   1.896 0.057898 .  
## age35-49             0.1435643  0.3101981   0.463 0.643497    
## age50-64            -0.2696730  0.2772296  -0.973 0.330681    
## age65-79            -1.0277636  0.2922478  -3.517 0.000437 ***
## age80+              -2.0710935  0.3598878  -5.755 8.67e-09 ***
## sexmasc             -0.2519979  0.1314381  -1.917 0.055208 .  
## income1             -0.0154454  0.0276949  -0.558 0.577051    
## hea1                -0.0720689  0.0064107 -11.242  < 2e-16 ***
## mar_1Yes            -0.5831916  0.3266439  -1.785 0.074196 .  
## score_lon1           0.3268351  0.0406748   8.035 9.33e-16 ***
## income1:hea1        -0.0001597  0.0005689  -0.281 0.778968    
## mar_1Yes:score_lon1  0.1290971  0.0668269   1.932 0.053383 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2585.4  on 3566  degrees of freedom
## Residual deviance: 1937.0  on 3555  degrees of freedom
##   (1186 observations deleted due to missingness)
## AIC: 1961
## 
## Number of Fisher Scoring iterations: 6

# Coefficients as odds-ratios
exp(coef(fit))

##         (Intercept)            age35-49            age50-64            age65-79              age80+             sexmasc             income1                hea1            mar_1Yes          score_lon1        income1:hea1 mar_1Yes:score_lon1 
##           2.4663996           1.1543810           0.7636291           0.3578063           0.1260479           0.7772463           0.9846733           0.9304668           0.5581142           1.3865729           0.9998403           1.1378006

plot(fit)

Gràficas de modelos lineales simples usando `ggplot`: una única variable predictora

fit <- lm(hea1 ~ income1, data = dataw1)
dataw1 %>%
  filter(!is.na(hea1), !is.na(income1)) %>% # nos quedamos con las observaciones que no tienen missings en el modelo
  # filter(!if_any(c(hea1, income1), is.na)) %>% # lo mismo de forma resumida
  cbind(predict(fit, interval = "confidence")) %>% # añadimos la predicción del modelo y los intervalos de confianza como nuevas columnas
  ggplot(aes(x = income1)) + # creamos el cuadro y añadimos el eje x
  geom_point(aes(y = hea1)) + # añadimos los punos d la base de datos
  geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .15) + # añadimos los CI; también es posible como sigue
  # geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
  # geom_line(aes(y = upr), color = "red", linetype = "dashed") +
  annotate("text", 
           label = paste("La correlación es", round(cor(dataw1$hea1, dataw1$income1, use = "complete.obs"), 2)),
           x = 30, y = 15) +
  theme_bw() # cambiamos el tema por defecto

Gràficas de modelos lineales usando `ggplot`: una variable predictora continua y otra categórica

library(ggplot2)
fit <- lm(hea1 ~ income1 + edu1, data = dataw1)
dataw1 %>%
  # filter(!is.na(hea1), !is.na(income1), !is.na(edu1)) %>% # nos quedamos con las observaciones que no tienen missings en el modelo
  filter(!if_any(c(hea1, income1, edu1), is.na)) %>% # lo mismo de forma resumida
  cbind(predict(fit, interval = "confidence")) %>% # añadimos la predicción del modelo y los intervalos de confianza como nuevas columnas
  ggplot(aes(x = income1, color = edu1)) + # creamos el cuadro y añadimos el eje x
  geom_point(aes(y = hea1)) + # añadimos los punos d la base de datos
  geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
  geom_ribbon(aes(ymin = lwr, ymax = upr, color = NULL, fill = edu1), alpha = .15) # añadimos los CI

Gràficas de modelos lineales generalizados simples usando `ggplot`: una única variable predictora

library(ggplot2)
fit <- glm(dep1 ~ income1, data = dataw1, family = binomial)
dataw1 %>%
  filter(!if_any(c(dep1, income1), is.na)) %>%  # nos quedamos con las observaciones que no tienen missings en el modelo
  mutate(fit = predict(fit, type = "response")) %>% # añadimos la predicción del modelo como nuevas columnas
  ggplot(aes(x = income1)) + # creamos el cuadro y añadimos el eje x
  geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
  labs(y = "Predicted probability of dep1 == Yes", 
       title = "Predicting  depression from income through a logistic model",
       caption = "Data source: A cohort study") + # especificamos la etiqueta del eje y, del título y de pie de figura
  theme_classic() #cambiamos el tema por defecto

Gràficas de modelos lineales generalizados usando `ggplot`: una variable predictora continua y otra categórica

library(ggplot2)
fit <- glm(dep1 ~ income1 + edu1, data = dataw1, family = binomial)
dataw1 %>%
  filter(!if_any(c(dep1, income1, edu1), is.na)) %>%  # nos quedamos con las observaciones que no tienen missings en el modelo
  mutate(fit = predict(fit, type = "response")) %>% # añadimos la predicción del modelo como nuevas columnas
  ggplot(aes(x = income1, color = edu1)) + # creamos el cuadro y añadimos el eje x
  geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
  labs(y = "Predicted probability of dep1 == Yes") + # especificamos la etiqueta del eje y
  theme_minimal() # cambiamos el tema por defecto

Otras formas de obtener resultados estadísticos en R; Ejercicios

De la misma manera que tidyverse hace referencia a una colección de librerías que simplifican la manipulación de datframes, existe otra colección de librerías, easystats (https://easystats.github.io/easystats/), cuyo objetivo es mostrar los resultados de tests y modelos estadísticos con mejor formato.

La librería modelsummary dispone de 2 funciones que permiten visualizar multiples resultados en tablas bien organizadas:

modelsummary (https://vincentarelbundock.github.io/modelsummary/articles/modelsummary.html) y
datasummary (https://vincentarelbundock.github.io/modelsummary/articles/datasummary.html).

Otras librerías que permiten visualizar descriptivos y resultados resumidos en tablas son:

summarytools (https://cran.r-project.org/web/packages/summarytools/vignettes/introduction.html) y
compareGroups (https://cran.r-project.org/web/packages/compareGroups/vignettes/compareGroups_vignette.html).

Existen distintas librerías que facilitan la representación gráfica de modelos, predicciones, …:

jtools (https://cran.r-project.org/web/packages/jtools/vignettes/effect_plot.html)
effects (https://cran.r-project.org/web/packages/effects/vignettes/predictor-effects-gallery.pdf)
ggeffects(https://strengejacke.github.io/ggeffects/)

Ejercicios:

Consultar la ayuda de cor. Calcular la correlación de Spearman. Cómo tendríamos que hacer si hay valores missing?
Instalar y cargar las librerías performance y parameters y probar las funciones model_performance y model_parameters aplicadas a fit.

Otras funciones y librerías útiles en R

A veces conviene guardar el nombre de un objeto en una variable y trabajar con la misma.

x <- 4       # asignamos el nombre x al objeto 4
w <- "x"     # guardamos el nombre x en la variable w
x

## [1] 4

## [1] "x"

get(w)       # con la función get recuperamos el objeto asignado al nombre (x) guardado en la variable w

## [1] 4

assign(w,5)  # asignamos el valor 5 al caracter guardado en w (x)
get(w)

## [1] 5

Si queremos ejecutar el código en R guardado en un archivo (habitualmente con la extensión .r), se puede ejecutar a través de la función source. Por ejemplo source("fichero.r").
En ocasiones nos puede interesar el número de caracteres en una cadena. Lo podemos obtener mediante la función nchar.

nchar(paste(x, 4)) == 3

## [1] TRUE

Tal como la librería haven guarda las etiquetas de los objetos como atributos de los mismos, a los que se puede acceder por medio de la función attributes otras librerías como labelled, sjmisc o sjlabelled también tienen múltiples funciones creadas para trabajar con etiquetas.
La función de R sample permite seleccionar elementos de un vector aleatoriamente. La función slice_sample de la librería dplyr permite seleccionar filas (observaciones) al azar de un data frame.
La función ggsave permite guardar gráficas generadas a partir de ggplot, especificando el formato, medida y calidad entre otros.

Apéndice

Ojo:
- Hay una cantidad enorme de librerías, muchas de las cuales comparten funciones (porque unas llaman a otras), pero también muchas de las cuales comparten nombres de funciones que son distintas (como muestra la advertencia al cargar la librería con funciones que ya están en el entorno). Por ejemplo:
- Las librerías plyr y dplyr comparten varias funciones, entre ellas rename. Si se carga primero dplyr y después plyr R asocia rename a la función de plyr, pero al usarla como la función de dplyr, es probable que se produzca un error.
- Es de los más frecuentes y extraños, por ello hay varias maneras de solucionarlo. Esencialmente, mediante el operador ::

library(dplyr); library(plyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

# ?rename # con la ayuda podemos ver que esta función está en varias librerías actualmente cargadas
names(ac1)

## [1] "q0002_hhid" "number_id"  "ID_ECS"     "sex"        "age"        "mar1"       "edu1"       "hea1"       "grups"

names(rename(ac1, hea = hea1)) # error

## Error in rename(ac1, hea = hea1): unused argument (hea = hea1)

names(dplyr::rename(ac1, hea = hea1)) # especificamos que la función que queremos es la de dplyr

## [1] "q0002_hhid" "number_id"  "ID_ECS"     "sex"        "age"        "mar1"       "edu1"       "hea"        "grups"

names(plyr::rename(ac1, hea = hea1)) # comprobamos que la función que estaba usando antes era la de plyr, ya que genera el mismo error

## Error in plyr::rename(ac1, hea = hea1): unused argument (hea = hea1)

rename <- dplyr::rename # otra solución es asignar al nombre la función de dplyr
names(rename(ac1, hea = hea1)) # ya no hay error

## [1] "q0002_hhid" "number_id"  "ID_ECS"     "sex"        "age"        "mar1"       "edu1"       "hea"        "grups"

Algunos recursos

Tutoriales introductorios a R:

Blogs:

https://statisticsglobe.com/r-programming-language/
enfocados a la estadística, pero no sólo (de hecho, diría que los siguientes son los mejores):
- http://www.sthda.com/english/
- http://www.flutterbys.com.au/stats/course.html
- http://www.r-tutor.com/ (incluye contenido de estadística bayesiana y machine learning, por otra parte un poco anticuado)
enfocado a la investigación en psicología (y vinculado a la librería psych):
- http://personality-project.org/r/
https://www.statmethods.net/index.html

Multitud de libros:

https://bookdown.org/
https://www.bigbookofr.com/ (contenidos parecidos a los del link anterior, en menor cantidad, aunque clasificados por materias)

Cursos online (MOOC’s) en plataformas como Coursera, edX, etc.

Más recursos

Acerca de librerías o materias específicas:

sobre librerías en general:
- https://rdrr.io/
- https://cran.r-project.org/web/views/
gráficas en general:
- https://www.r-graph-gallery.com/
- https://r-charts.com
- https://ggplot2-book.org, el libro de ggplot2
- https://exts.ggplot2.tidyverse.org/gallery/, librerías que extienden las funcionalidades de ggplot2
formas de mostrar descriptivos:
- https://dabblingwithdata.wordpress.com/2018/01/02/my-favourite-r-package-for-summarising-data/
meta-análisis:
- https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/
neuroimagen:

Recursos de RStudio

El foro https://community.rstudio.com/ principalmente centrado (aunque no exclusivamente) en cuestiones relacionadas con RStudio y con las librerías de la familia tidyverse
La página web https://support.rstudio.com/hc/en-us contiene material acerca de RStudio
El blog https://blog.rstudio.com publica contenido diverso relacionado con los productos de la compañía RStudio
El blog https://rviews.rstudio.com está más enfocado al usuario de R, aportando información y ejemplos del uso de distintas librerías, así como acerca de conferencias de usuarios de R
El blog https://www.tidyverse.org/blog/ publica habitualmente entradas acerca de las nuevas librerías de la familia tidyverse o de actualizaciones relevantes

Podéis consultarme vustras dudas a través de mi dirección de email: iago.gine@sjd.es

Introducción a R

The R Project for Statistical Computing

R como calculadora

Necesitamos un editor

Algunas nociones sobre RStudio

Qué es lo primero que nos encontramos en R? Objetos, funciones, comentarios, (directorios) y librerías

Siguientes pasos: formas de obtener ayuda, argumentos de una función y abrir un fichero con formato de Excel

Abrir bases de datos y ejecutar funciones: R no es como Stata.

Objetos complejos. Clases. Bases de datos

Otras características y algunas funciones prácticas

Primeros pasos con una base de datos.

Acceso a las variables de una base de datos: los operadores $ y [[

Otras funciones y operadores lógicos y aritméticos en R

Ejercicios

Transformación de variables. Factores

Transformación de variables. Factores ordenados

Transformación de una variable en una base de datos (con la librería dplyr)

Pipe

Transformación de varias variables en una base de datos con la librería dplyr

Otras operaciones con funciones de la librería dplyr: selección de variables

Otras operaciones con funciones de la librería dplyr y funciones de R para el cálculo de estadísticos

Cálculo de prevalencias:

Vemos las variables y guardamos los resultados

Ejercicios

Condiciones y loops

Funciones *apply

Variables y dimensiones de las bbdd

Fusión (merge) de bbdd

Bases de datos en formato de Stata

Variables con etiquetas

Combinar bases de datos

Otras opciones para combinar bases de datos

Rellenando datos

Ejercicios

Pivotar (vertical/largo –> horizontal/ancho)

Pivotar (horizontal/ancho –> vertical/largo)

Correlación y tests de normalidad

Tests chi-cuadrado y de Fisher

Tests univariantes

Prueba t de Student (o t-test)

Prueba Wilcoxon-Mann-Whitney (también test U de Mann-Whitney o test de la suma de rangos de Wilcoxon)

ANOVA

Regresión lineal

Regresión logística

Gràficas de modelos lineales simples usando ggplot: una única variable predictora

Gràficas de modelos lineales usando ggplot: una variable predictora continua y otra categórica

Gràficas de modelos lineales generalizados simples usando ggplot: una única variable predictora

Gràficas de modelos lineales generalizados usando ggplot: una variable predictora continua y otra categórica

Otras formas de obtener resultados estadísticos en R; Ejercicios

Otras funciones y librerías útiles en R

Apéndice

Algunos recursos

Más recursos

Acceso a las variables de una base de datos: los operadores `$` y `[[`

Transformación de una variable en una base de datos (con la librería `dplyr`)

Transformación de varias variables en una base de datos con la librería `dplyr`

Otras operaciones con funciones de la librería `dplyr`: selección de variables

Otras operaciones con funciones de la librería `dplyr` y funciones de R para el cálculo de estadísticos

Gràficas de modelos lineales simples usando `ggplot`: una única variable predictora

Gràficas de modelos lineales usando `ggplot`: una variable predictora continua y otra categórica

Gràficas de modelos lineales generalizados simples usando `ggplot`: una única variable predictora

Gràficas de modelos lineales generalizados usando `ggplot`: una variable predictora continua y otra categórica